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Abstract 

A time- accurate, upwind, finite volume method for computing compressible flows on unstructured 
grids is presented. The method is second order accurate in space and time and yields high resolution 
in the presence of discontinuities. For efficiency, the Roe approximate Riemann solver with an entropy 
correction is employed. In the basic Euler/Navier- Stokes scheme, many concepts of high order upwind 
schemes are adopted: the surface flux integrals are carefully treated, a Cauchy-Kowalewski time- stepping 
scheme is used in the time-marching stage, and a multidimensional limiter is applied in the reconstruction 
stage. 

However, even with these up-to-date improvements, and a Roe Riemann solver with entropy correc- 
tion, the basic upwind scheme is still plagued by the so-called “pathological behaviors,” e.g., the carbuncle 
phenomenon, the expansion shock, etc. A systematic solution to these limitations is presented, which uses 
a simple dissipation model while still preserving second order accuracy. The modified, stabilized scheme 
is referred to as the enhanced time-accurate upwind (ETAU) scheme in this paper. 

The unstructured grid capability renders flexibility for use in complex geometry; and the present ETAU 
Euler/Navier- Stokes scheme is capable of handling a broad spectrum of flow regimes from high supersonic 
to subsonic at very low Mach number, appropriate for both CFD (computational fluid dynamics) and CAA 
(computational aeroacoustics). Numerous examples are included to demonstrate the robustness of the 
method. 

1 Introduction 

Finite volume (FV) schemes are gaining popularity in computational fluid dynamics (CFD) and computational 
aeroacoustics (CAA) primarily due to their robustness and geometric flexibility. The FV schemes are based on the 
Gauss divergence theorem applied to a control volume (CV) and consist of two steps. In the first (reconstruction) step, 
with given initial conditions, the cell average flow variables are reconstructed into linear or higher order polynomials 
within the CV. The second (evolution) step involves computing the surface fluxes of the CV; the cell averaged values of 
flow variables are then obtained for a solution at the next time level. The surface flux calculation in these schemes can 
be categorized into two types: the centered schemes and the upwind schemes [1]. While the upwind schemes require 
a Riemann solver (exact or approximate), the centered schemes, such as the NT (Nassayahu-Tadmor) [2] scheme and 
the CE/SE scheme [3, 4], do not. 
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There are currently many popular upwind schemes being used in CFD and CAA. The Godunov scheme and the 
TVD high resolution schemes [5, 6] are considered to be the fundamental upwind schemes. In the past decades there 
has been significant progress in improving the accuracy of the upwind methods by using higher order approximations, 
e.g., the ENO (essentially non-oscillatory) [7] and WENO (weighted ENO) [8] schemes; the DG (discontinuous 
Galerkin) [9] scheme; the SV (spectral finite volume) and the SD (spectral difference) schemes [10, 11]. 

The purpose of this paper is to present a practical method that avoids high computer costs but with reasonable 
accuracy while retaining its robustness. In this work a basic upwind scheme is chosen similar to the one presented, 
for solving the Euler equations, in [1]. As the upwind methods are prone to exhibit “pathological behaviors” [12- 16] 
such as the carbuncle phenomenon and the expansion shock, the method presented here will provide a cure for these 
undesired phenomena. 

In the present paper, we focus on the construction of an enhanced time-accurate upwind (ETAU) FV scheme for 
the Euler/Navier- Stokes equations on unstructured grids (triangular or tetrahedral), which removes the pathological 
behaviors through a simple dissipation model. Other researchers have considered the role and implementation of 
dissipation models in upwind schemes [1, 12 - 16]. However, the present work takes a new systematic approach to 
numerical dissipation to broaden the capabilities of the basic upwind scheme. The ETAU scheme is second order 
accurate in space and time, and stable for long run times. It is capable of computing flows over a wide range of flow 
regimes which include discontinuities as well as very low Mach number viscous flows, making the scheme appropriate 
for practical CFD and CAA problems. 

The governing equations and the basic unstructured Euler/Navier- Stokes solver are briefly described in Section 2, 
with simple implementations of the boundary conditions, in particular, the non-reflecting boundary condition (NRBC). 
Section 3 is devoted to the cure of the pathological behaviors that warrants a stable scheme - ETAU. Several examples 
are presented in Section 4 that demonstrate how such annoying behaviors are overcome in a systematic way. In Section 
5, the ETAU Euler/Navier- Stokes solver is tested in numerical examples from supersonic to low subsonic flow speed, 
with emphasis on stable, long run time viscous flows for CAA computations. Concluding remarks are addressed in 
Section 6. 

2 The Basic Unstructured Euler/Navier-Stokes FV Solver 

Because computer time and memory still need to be considered when computing time-accurate unsteady flows, 
an upwind FV scheme that is second order accurate (in the absence of flow discontinuities) in space and time is 
chosen. The basic upwind scheme includes many attributes often found in the construction of high order schemes: 
the geometries are calculated as accurately as possible, the use of a Cauchy-Kowalewski type evaluation of the time 
derivative, and a multidimensional limiter for the reconstruction stage. A brief description of the basic upwind scheme 
is outlined below, beginning with the Euler/Navier-Stokes (E/N-S) equations in conservation form. 

2.1 Conservation form of the unsteady compressible Euler/Navier-Stokes equations 

The conservation form of the two-dimensional E/N-S equations are briefly sketched here. The three-dimensional 
E/N-S equations can be treated in a similar way. Let p, u , v, p , and y be the density, the two velocity components, static 
pressure, and constant specific heat ratio, respectively. The two-dimensional unsteady E/N-S equations are written in 
the standard conservation form: 

U,+F x + G, = 0, (1) 

where x, y, and t are the streamwise and transverse coordinates and time, respectively. The flux vectors are further 
split into inviscid and viscous fluxes: 

F = Fj — F v , G = G* — G v . 

The conservative flow variable vector U, and the inviscid flux vectors F ? and G ? are given in nondimensional form as: 
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Figure 1 . (a) A typical unstructured triangular grid in two-dimensional space, (b) control volume in £ 3 , (c) ghost cell is a mirror image of the 

boundary cell. 


Here the internal energy has the form e = + 1/2 (w 2 + v 2 ), and the enthalpy is H = p/p + e. The nondimension- 

alized viscous flux vectors F v and G v are written as: 

/ 0 \ 

_ n(2u x — §V-V) 

V K V X + Uy) ’ 

+ (i Uy + v*)v - l (V-V )u+ Tr^Th~ )] / 

and 

/ 0 \ 
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where u x ,u y ,v x , are the flow velocity gradients with Pr being the Prandtl number, /n = l /Re is the dynamic viscosity 
where Re is the Reynolds number; the velocity divergence is V-V = u x + v y . For air at standard conditions, Pr = 0.72 
and y = 1.4. 


2.2 The time-marching (evolution) stage 

A typical two-dimensional unstructured triangular grid cell used in the scheme is illustrated in Fig. 1 . Here, A ABC 
is a triangular cell centered at O and D,E,F are the centers of the neighboring triangular cells (compact stencil). The 
flow variables at the previous time step are stored at these triangle cell centers. By considering (x,y, t) as coordinates of 
a three-dimensional Euclidean space, £ 3 , and using the Gauss divergence theorem, it follows that Eq. (1) is equivalent 
to the integral conservation law: 

^I m *ds = 0, m=l,2,3,4, (2) 

where S denotes the surface around a space-time control volume (CV) in £3 (e.g. the prism ABC — A'B'C' in Fig. lb) 
and I m = (F m , G m , U m ). For the surface integrals in Eq. (2), according to the philosophy of high order upwind schemes 
[10, 11, 17,], gaussian quadratures are required. Here, for second order accuracy, the integrands are approximated by 
linear functions and the only gaussian point is the centroid of each space-time hyper-surface (e.g., M in Fig. lb, or in 
Fig. 2a for a tetrahedral cell). The simple quadrature of the surface integrals in Eq. (2) leads to the update formulation 
of U at O from time step n to n + 1 : 

U«+1 = U ■" "A £ [F l +l/2 (n x ) k + of 1 ,/2 (ny)t] Al k , (3) 

k= 1 

where F&, G&, (n x )k, ( n y )k , k= 1,2,3, are respectively the two flux vectors, the out-going unit normal vector com- 
ponents at the centers of the three cylindrical surfaces ABB' A ' , BCC'B' and CAA'C ' . The lengths of the edges of the 
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Figure 2. (a) A typical unstructured tetrahedral grid cell ABCD in three-dimensional space, with O as its centroid, (b) ghost cell is a mirror image 

of the boundary cell, centered atP. 


triangular cell A ABC, AB , BC and CA, are represented by A 4, k— 1,2,3; As the area of A ABC and At the time step 
size. 

The inviscid part of the fluxes F^ +1 ^ 2 , Q^ +1 / 2 are generated via the Riemann solver. Here, the Roe approximate 
Riemann solver [18] with an entropy correction [1] is employed and works very well. At time level t = nAt , the 
conservative flow variable vector U n is given at all triangular cell centers (e.g. 0,D,E,F in Fig. la or P in Fig. 2a for 
three-dimensional space). The spatial gradients U x and can then be reconstructed via a multidimensional limiter 
(as described below), and conveniently stored for further use. The left and right (L and R in Fig. 1 and Fig. 2) states 
at the center of each surface of the space-time CV are established by linear Taylor expansion from their corresponding 
cell centers. For instance, at the center, M, of the surface ABB' A' (Fig. lb, Fig. 2a), the R and L states are respectively: 

uf 1/2 = V n D + (U”) D Ax+ (U>Ay+ (U f >y , (4) 

U" +1/2 = U n 0 + (V n x ) 0 Ax+ (U>Ay + (U f >f • (5) 

Here, the subscripts O and D indicate which cell center flow variables are used, Ax and Ay correspond to either DM or 
OM. Note that since the out-going unit normal is chosen, the L-state is always evaluated from the current cell center 
O and the R-state always from its neighboring cell centers (£>,£, or F). 

The procedure that is outlined by Huynh [1] for solving the Euler equations is used in this work for the solution of 
the inviscid part of the Navier- Stokes equations. Huynh also showed that the Euler scheme is second order accurate in 
space and time. This is used as the basis of a flow solver that will be augmented later to cure the undesired pathological 
behaviors of this basic upwind scheme. 

The evaluation of the viscous fluxes F v and G v follows the formulation in Section 2.1, but only the R-state data is 
used. As the flow is assumed to be linear within each cell, spatial derivatives are constant therein and can be applied 
directly without Taylor expansion. 

In the time direction, the time derivative U* is evaluated at time step nAt as described below. 


2.3 The Cauchy-Kowalewski/Lax-Wendroff time stepping 

The Cauchy-Kowalewski concept makes use of the governing equations to represent the time derivative by the 
spatial derivatives. Lax and Wendroff were the first to apply the concept to the Euler equations in their well-known 
Lax-Wendroff scheme. The idea was followed by other researchers, e.g., the NT and CE/SE schemes [2,4] use this 
technique to obtain U” at t — nAt. Recently, Toro further extended the idea in deriving a solution for the GRP (gen- 
eralized Riemann problem) in their ADER scheme [19]. Here, for second order accuracy in time, evaluation of U* 
follows the Lax-Wendroff procedure: 


U« = — F" 


r*n 


( 6 ) 


As the Reynolds number, Re, is usually high (or viscosity, p, small), F and G can be well-approximated by their 
inviscid portions, (6) is simplified to: 


U n = -(dfiyu" _ (dfh) 
' 1 au } * 1 au ’ 


! U" 


(7) 
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U? is then substituted respectively to (4) or (5) to obtain the initial right and the left states U^ +1 ^ 2 and U ^ +1//2 . Primitive 

variable equivalents of U^ +1//2 and X} n ^ 1 ^ are then provided as R and L states to the Roe approximate Riemann solver 
[1, 18]. A numerical flux for a given face is thus evaluated. Once the flux data is computed at each face, the averaged 
solution is updated to the (n + 1)A t time level (the top face of the CV, A A'B'C') using equation (3). 


2.4 Evaluation of spatial gradients by a multidimensional limiter and the reconstruction stage 

The reconstruction stage uses averaged flow data U updated at the current cell center O and at its three neighboring 
cell centers D,E, and F (Fig. la). In general, the reconstruction process requires solving an overdetermined system 
(see, e.g.,[17]). For second order accuracy, U must be reconstructed into a linear vector function within the current 
cell via a limiter. Here, a multidimensional as opposed to dimension by dimension limiter is employed. Let u denote 
a component of U and / = \OD\; by simple finite difference, a linear equation for the gradients ( u x , u y ) is obtained: 


ud — uo 

7 


_ xd-xo . yD-yo 

i Uy i 


or 


(xd - Xo)u x + (yD - yo)u y =U D -Uo- 


Here, the subscripts O and D denote the corresponding cell centers. Similarly, two more such linear equations are 
available from the u data at E and O as well as at F and O. 


(x E -xo)u x + (y E —yo)u y = U E - uo , 


(x F -xo)u x + (y F - yo)u y = u F - u 0 - 

Generally, any combination of two of these equations yields a set of spatial gradients for u. There are three sets of 

such gradients, namely, (u^ x \u^), (u^ x \u^) and (u^ x \u^). Their moduli or I 2 norms, m; = \J ( u x ^) 2 + (w^) 2 , 
i = 1,2,3, are calculated as a measure. Then, a multidimensional limiter based on the magnitude of the norms is 
applied to all spatial gradients to achieve a single, unified set of gradients for U. Two such multidimensional limiters 
are recommended here: the minmod limiter and the extended van Albada limiter [20] (weighted averaging). 

When the minmod limiter is employed, the gradient {u x \uf > ) corresponding to m ? - = mm (mi, m 2 , m 3 ) is chosen. 
When the extended van Albada limiter is used, let w\ = ,w 2 = (mim 3 ) a ,W 3 = (mim 2 ) a , a > 0, the final 

gradients are evaluated through weighted averaging: 

( 1 ) , ( 2 ) . ( 3 ) ( 1 ) . ( 2 ) . ( 3 ) 

_ W\Ux +W2Ux J +W3Ux _ _WiUy J +W2Uy J +W3Uy' 

U X , Uy . 

W\ + W2 + W 3 W\ + W2 + W 3 

The extended van Albada limiter or weighted averaging was previously used in the CE/SE method [3, 4]. Different 
as provide numerical dissipation at different levels. A small a, e.g.,a = 0.5, usually yields less dissipation, and is 
appropriate for aeroacoustics computations. In the presence of a shock or a contact discontinuity, a larger a(= 2.0) is 
needed. In the rare extreme situation with high Mach number and strong shocks, the computed gradients u x , u y need 
to be further limited by a factor in order to keep the numerical procedure stable. 

After the gradients for each component of U are evaluated, the reconstructed U is linear within the entire grid cell, 
A ABC (Fig. la and b). The solution can be considered as either second order accurate or high resolution. An impor- 
tant advantage of the multidimensional reconstruction is that a simple and robust absorbing nonreflecting boundary 
condition (NRBC) can be applied, as sketched in the next subsection. 

2.5 Boundary conditions 

The solver is completed with treatments of some frequently used boundary conditions. For convenience, all the 
ghost cells are generated as a mirror image of their corresponding boundary cell (A ABC in Fig. lc). The given physical 
boundary conditions (BCs) are specified at the ghost cell centers (e.g. D in Fig. lc). 
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2.5.1 slip wall BCs For a slip wall boundary, flow variables U = (C/i,C/2, U^U^) 7 at the ghost cell center D 
(Fig. lc) are mirror images of those at the boundary cell center O. Let W = (C/ 2 , C/ 3 ), 

( C/i )d = ( C / i ) o , W/) =Wo-2(Wo*n)n, (Ua)d = (Ua)oi 

where n is the out-going unit normal at the boundary (A# in Fig. lc). Spatial gradients at the ghost cell center also 
need to be defined in an appropriate way. 

2.5.2 no-slip wall BCs For a no-slip wall, the BC is slightly different. 

(Ui) D = (Ui)o, W d = -W o, (C/ 4 )d = ( 1 / 4 ) 0 . 

2.5.3 nonreflecting boundary conditions (NRBCs) One of the most important advantages of the mul- 
tidimensional reconstruction is that it allows a simple but robust absorbing NRBC. In principle, the NRBC is based on 
the decomposition of the local Euler solution into its Fourier modes - plane waves. A C 1 continuity criterion for the 
flow variables is inferred as a generic NRBC [21, 22]. After the multidimensional reconstruction, U is linear and hence 
C 1 continuous within the CV , in particular, across the boundary surface ( dash line AB in Fig. lc). Then practically 
no reflection is generated, details can be found in [22]. The NRBC can be consistently extended to three-dimensional 
space, only the line AB now becomes a boundary surface element A ABC (Fig. 2b). In principle, the difference between 
the conventional characteristic-based NRBC and the present NRBC is, the former is based on decomposition of the 
Euler equations into characteristic fields, and the latter is based on a local decomposition of the solution into plane 
waves or Fourier modes. But in practice, no such decomposition operation is required. Any given physical boundary 
conditions at the ghost cell center (e.g. D in Fig. 2a) are automatically nonreflecting , no additional numerical im- 
plementation is needed. Often, when discrepancy between the BC and the flow in the domain interior is developed, 
an absorbing layer is consequently formed around the boundary as a buffer zone. The absorbing layer is somewhat 
similar to the matched layer in the PML (perfectly matched layer) method [23] for NRBC, but they are two different 
methods. More details on the NRBCs can be found in [21, 22]. 


2.6 Extension to three-dimensional space 

The basic E/N-S upwind scheme can be consistently extended to three-dimensional space. The governing equa- 
tion, Eq.(l), is now extended to 

U,+F x + G, + H z = 0. (8) 

Here, similar to Section 2.1, U, F, G, and H are the three-dimensional versions of the conservative flow variables and 
the flux vectors. All the numerical implementations from Sections 2. 2-2. 5 can be extended to the three-dimensional 
case. The update formulation, Eq. (3), becomes 


V 


n + 1 _ 1 


= U" - f v £ K +1/2 (n x ) k + G n k +l/2 (n y ) k + K n k +1/2 (n z ) k ]As k , 


k= 1 


( 9 ) 


where Av and A Sk, k= 1,2, 3,4, represent the volume and the surface areas of the current tetrahedral cell. Figure 2a 
shows the tetrahedral cell ABCD centered at O and one of its four neighbors centered at P. For a boundary cell ABCD , 
the corresponding ghost cell is its mirror image (Fig. 2b). 


3 Cure of the Pathological Behaviors 

The basic scheme described above is an improved upwind scheme of the Godunov type. The scheme works well 
with many benchmark type problems but still suffers from the general symptom of the so-called pathological behaviors 
described by Quirk [12] and others. The pathological behaviors and their cure are discussed below. 

3.1 The pathological behaviors and how they occur 

Despite the great success of the upwind schemes, there are still some cases of failure reported by researchers. For 
example, the carbuncle phenomenon in the bow shock of a blunt body, the expansion shock when a supersonic flow of 
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Figure 3. (a) and (b): A rectangular cell ABCD near a bow shock of a blunt body; (a) the cell ABCD, showing the forming of dissipation in a 

Riemann solver, L and R correspond to left and right states; (b) location of the cell relative to the shock ; (c) grid structure near a sharp edge in the 
expansion shock case, P and Q are cell centers . 


high Mach number passes over a sharp edge, and the kinked Mach stem when a supersonic flow passes a wedged wall, 
etc. [13 - 16]. These pathological behaviors occur mostly when a grid line is aligned with a strong shock or expansion. 
It is generally agreed that insufficient dissipation and the consequent local numerical instability lead to such failures 
[12, 13, 15]. 

In the literature, for the Godunov type upwind schemes, a general cure is switching to different Riemann solvers 
to gain some more dissipation. For example, an expansion shock can be cured by switching from a Roe approximate 
Riemann solver to an HLLE Riemann solver [12]. A comparison of the performances of various Riemann solvers or 
flux- splitting is given in [13]. One needs to know in advance which Riemann solver is appropriate for the particular 
flow problem being considered. This reduces the robustness of the scheme. Lin [14] provides a way to add dissipation 
to FDS (flux difference splitting) schemes to treat the carbuncle phenomena and the slow moving shock instability. 
Other curing measures may be found in the references of these papers. The recent residue distribution (RD) schemes 
(e.g. [24]) utilize a pointwise representation of the solution. As in a finite difference scheme, the unknowns are updated 
at the cell vertices and the stabilization mechanism is similar to artificial viscosity. 

In his investigation, Xu [15, 16] further argued that for a Godunov type upwind scheme, 

The dissipation required by the numerical stability mainly comes from the Riemann solver at the cell interface. 
With the given L and R states across the cell interface (L and R in Fig. 3a), the Riemann solver attempts to 
compromise the two different states to an intermediate state for the flux by some kind of averaging. For example, 
in the Roe approximate Riemann solver, Pl/(Pl + ^/PlPr) is used as the weighing factor in the weighted 
averaging for the flow variables. The averaging process is accompanied by an entropy increase and numerical 
dissipation is generated. 

The amount of dissipation varies subject to the L and R states and is in the direction normal to the cell interface. 
The more they differ from each other, the more the numerical dissipation. 

After applying a Riemann solver to the interface between two adjacent cells, the divergence theorem stipulates 
that only the normal flux vector components at the interface are considered in the evaluation of surface fluxes 
and all the tangential components are ignored. In the direction tangent to the surface, no wave is assumed to 
occur and the numerical dissipation is absent. 

Based on Xu’s arguments, we are able to conduct some qualitative analysis to the pathological behaviors such as the 
carbuncle phenomena. For convenience, we assume a rectangular grid as shown in Fig. 3a and b. We also assume that 
the grid lines AD and BC are aligned with the bow shock. For the grid cell ABCD located near the stagnation point of 
the blunt body (Fig. 3a), as the flow states jump sharply across BC due to the bow shock (Fig. 3b), the Riemann solver 
(RS) at BC introduces numerical dissipation. The dissipation is directional, no dissipation in the tangential direction 
of BC occurs. For the surface (edge) CD , as the L and R states are almost identical (e.g., vp = vr), the RS provides 
little or no dissipation in its normal direction. Similar analysis can be applied to the other two surfaces AB and AD. 
Thus, the numerical solution is vulnerable to the local temporal or spatial instabilities in the direction tangent to the 
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Figure 4. Adding numerical dissipation to a one-dimensional scheme, (a) grid for an unconditionally unstable scheme; (b) grid for the Lax- 
Friedrichs scheme. 



Figure 5. Dissipation models in two-dimensional and three-dimensional spaces: (a) for two-dimensional triangular cell, O is the cell center, 
M.N.P are the cell edge centers; (b) for three-dimensional tetrahedral cell, O is the cell center, M.N.P. Q are the cell surface centers. 


shock, odd-even decoupling or (numerical) shear layer instability may take place and cause the carbuncle phenomena. 
A remedy is to add some extra numerical dissipation in the tangential direction to stabilize the numerical computation. 

Another typical pathological behavior is the expansion shock that forms in supersonic flow of high Mach number 
past a sharp edge (Fig. 3c). The expansion shock represents a different type of insufficient dissipation. Assume the 
grid line ABC is aligned with the edge body surface. Although the RSs at AB and BC generate numerical dissipation 
in the normal direction, the dissipation may not be enough, and an “expansion shock” still appears. Quirk [12] has 
suggested using a more dissipative HLLE Riemann solver and smearing the “expansion shock” to an expansion fan. 
Here again, the remedy is to add more numerical dissipation, but this time, in the direction normal to the cell interface. 

Analysis of the above two cases hints that if there exists a dissipation model which contributes additional dissi- 
pation in both tangential and normal directions or is omnidirectional, these pathological phenomenon will be cured. 
Details about this dissipation model are discussed below. 

For the unstructured grid we considered in the present paper, the symptom of local numerical instability may look 
less severe at first glance as the cell surface orientation is random. However, a similar situation may still occur along 
the boundaries. For example, it was found with the basic Euler/Navier- Stokes scheme, the no-slip wall condition does 
not work at all. The symptom was overcome after a dissipation model was used. 

For robustness and efficiency, we continue to use the Roe approximate Riemann solver (RS) in the present basic 
upwind scheme, and propose a simple dissipation model that is external to the RS. The model adds omnidirectional 
dissipation to and is embedded in the basic scheme. 


3.2 A sample one-dimensional dissipation model 

In searching for an appropriate dissipation model, we begin with a simple one-dimensional scalar advection 
equation 


du 

dt 



( 10 ) 
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Two simple one-dimensional schemes are reviewed. The first scheme is a scheme of forward difference in time and 
central difference in space (Fig. 4a) 

u n j+ l = w” + 0.5r(w^ +1 — u n j_ i) = u a , 

where r = At /Ax. It is well known that such a scheme is unconditional unstable. However, if u " is replaced by an 
average of u at the adjacent points, the scheme becomes the Lax-Friedrichs scheme (Fig. 4b), 

m" + 1 =0.5(w" +1 +M"_ 1 ) + 0.5r(M" +1 -M"_ 1 ) = u b , 

which is stable for r < 1. Here, we have learned that when replacing uj by an average of u at the adjacent nodes, a 
certain amount of dissipation is added to the scheme, turning the scheme from unstable to a stable one. Furthermore, 
in order to control the dissipation, Huynh [1] (ref. Eq. 3.11 therein) blends the two schemes by a weighted average: 

w" +1 =pM fl + (l-p)wfo, 

where 0 < P < 1 . This idea will be extended and applied to the multidimensional schemes below. 

3.3 Multidimensional artificial dissipation model 

It is straightforward to construct an omnidirectional dissipation model. For example, replacing XJ n in Eq.(3) by a 
weighted average of \J n at its neighboring cell centers D , E , and F as in Fig. la: 

U" = aiU£ + a 2 U£ + a 3 U£ 

will produce dissipation, where ai >0, a 2 > 0, a 3 > 0, and ai + a 2 + a 3 = 1. However, in order to control the 
amount of dissipation, a multidimensional dissipation model is proposed based on the above one-dimensional model. 
The process is fully embedded in the numerical procedure at no additional cost in operation. 

As shown in Fig. 5a, let A ABC be the current cell, M, N,P be respectively the mid-points of its edges AB,BC , and 
CA (i.e., cell surface centers). The triangle A MNP is also centered at O. The R states at these mid-points at time level 
n are extrapolated from their respective neighboring cell centers D,E,F: 

V n M = \J n D + (V n x ) D Ax + (V n y ) D Ay, 

U n N = Ul + (V'i) E Ax+(V n y ) E Ay, 

V n P mV n F + (V n x ) F Ax+(V n y ) F Ay. 

They are an intermediate step of Eq. (4) and cost no extra CPU time in computation. The average of these pointwise 
values 

U=(U^ + U£, + U£)/3 

is an approximation of U o with spatial smearing. Then, in (3) is replaced by a weighted average of \J n and U, i.e., 
Eq. (3) becomes 

U" +1 =(3U+(l-p)U"-^^[Ff 1/2 (^ + Gf 1/2 K) i ]A/ 4 , (11) 

^ k= 1 

Here (3, 0 < (3 < 1, is the weighing factor. By this weighted averaging process, some numerical dissipation is introduced 
to the scheme. Selection of the (3 values depends on the grid geometry and the flow. For most cases, the flow is 
insensitive to the value of (3 and (3 = 10 -4 — 10 -3 is appropriate for suppressing the pathological behavior. For flows 
with extremely high Mach number, P may need to be increased to 0. 1 — 0.2 or higher. Generally, setting P to 10 -3 is a 
good compromise. With the dissipation model added to the basic scheme, it is referred to the enhanced time-accurate 
upwind (ETAU) scheme. 

The ETAU scheme (11) falls within the class of second order accurate upwind schemes analyzed by Huynh [1] 
where Fourier analysis is conducted for weighted schemes on triangulated structured meshes. Further discussion on 
the accuracy of upwind finite volume schemes using arbitrary unstructured polyhedral grids for three-dimensional 
flows can be found in Delanaye and Liu [17]. 

Equation (11) can be consistently extended to a three-dimensional scheme with tetrahedral grids as shown in 
Fig. 5b. In the following, the cure of pathological behaviors is demonstrated in various examples. 
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Figure 6. Cure of expansion shock for a supersonic flow diffraction problem. (1) Godunov scheme with no cure (p = 0), the expansion shock is 
clearly shown; (2) Godunov scheme with cure (ETAU, (3 = 0.001), the expansion shock disappears; (3) high-resolution with minmod limiter; (4) 
high resolution with van Albada limiter. 


4 Examples for Curing the Pathological Behaviors 

In this section, we demonstrate how several well-known pathological behaviors are overcome with the present 
dissipation model. The general performance of the ETAU scheme will be illustrated in the next section. 

4.1 The expansion shock 

Consider the problem of a strong shock diffracting over a 90° edge [12]. Here, the shock Mach number is 
M s = 5.09. There are about 14, 400 triangular cells in the domain. In order that the expansion could occur, a horizontal 
grid line is imposed by extending the top solid wall surface line into the domain. Initially, the flow is set to the quiescent 
ambient condition: 

(po,wo,vo,po) = (1,0,0, 1/y), Y= 1.4, 
and M s = 5.09 conditions are imposed at the inlet: 

(p i,u u vuPi) = (5.0294,4.0779,0,21.4710). 

At the top, bottom, and along the surfaces of the rectangular block, the slip wall condition is imposed. A simple 
extrapolation condition is applied to the outlet boundary. Figure 6 shows the density contours with or without the cure. 
With (3 changing slightly from 0 to 0.001, the expansion shock disappears. 

4.2 The kinked Mach stem 

Another example of the pathological behaviors is the kinked Mach stem over a ramp of 30° [12]. The grid, though 
unstructured, is carefully designed with grid lines parallel to the ramp surface or to the y-axis to ensure occurrence of 
a kinked Mach stem. There are 36,240 triangular cells in the grid. Initially, the flow is set at the quiescent ambient 
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Figure 7. Supersonic flow over a ramp with recovery of the correct Mach stem: (1) Godunov scheme with no cure ((3 = 0), noise at the top 
wall and the kinked Mach stem are clearly shown; (2) Godunov scheme with cure (p = 0.001), the noise disappears and a correct Mach stem is 
recovered; (3) high-resolution with minmod limiter (p = 0.001 ); (4) high resolution with van Albada limiter (p = 0.001 ). 


condition: 

(po,Mo,vo,po) = (1,0,0, 1/y), 
and M s = 5.5 conditions are imposed at the inlet: 

(p i,Ui,v u pi) = (5.1489,4.4318,0,25.0893). 

Boundary conditions at the top, bottom, and the ramp surface are the slip- wall conditions. Simple extrapolation is 
imposed at the outlet boundary. Figure 7 demonstrates how a correct Mach stem is recovered by setting (3 = 10 -3 . 

4.3 The carbuncle phenomenon 

The carbuncle phenomenon occurs with most upwind schemes when computing supersonic flow past a blunt body 
with a bow shock. In this example, the blunt body is a half sphere (circle) with a slip-wall boundary condition. Simple 
extrapolation is applied to the outflow boundaries. There are about 80, 000 triangular cells in the computational domain 
and the grid lines are somewhat aligned with the bow shock to trigger the carbuncle phenomenon. 

As shown in Fig. 8, for a Godunov scheme (Plot a), the basic scheme ((3 = 0) with van Albada limiter (Plot c) 
or with a minmod limiter (Plot e), the carbuncle phenomenon still exists, even though an entropy correction has been 
imposed in the Roe approximate Riemann solver. However, the carbuncle phenomenon disappears when (3 = 0.01 is 
applied. Due the high Mach number, M = 10, the spatial gradients after limiting are further reduced by a factor of 0.1 
or 0.2 in order to maintain numerical stability. 

4.4 The slowly moving shock 

Here the one-dimensional problem is computed in a two-dimensional domain. A grid consisting of 800 uniform 
triangular cells is used, spanning between 0 < x < 32. Initially, a strong shock is located at x = 15 with the left states 
(p/,M/,/?/) = (3.86,0.81,10.34) and the right states (p r ,w r ,p r ) = (1,3.44,1). These states are also imposed as the 
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Figure 8. Supersonic M = 10 flow past a circular blunt body, showing carbuncle phenomena and their cure; (a) Godunov (1st order) scheme, 
(3 = 0, (b) Godunov scheme with (3 = 0.01 imposed; (c) upwind scheme with van Albada limiter, (3 = 0; (d) with (3 = 0.01 imposed; (e) upwind 
scheme with minmod limiter; (f) with P = 0.01 imposed. 



Figure 9. Suppression of spurious oscillation in one-dimensional slowly moving shock at t = 20 (10,000 steps); left: comparison for p = 0 and 
P = 0.3; right: comparison for p = 0.3 and p = 0.5, shows growing numerical dissipation with increased p. 


inflow and outflow boundary conditions. At the top and bottom of the domain, the usual reflective slip wall conditions 
are imposed. With At = 0.002, after running 10, 000 steps, the slow shock only moves from v = 15 to v = 17, as shown 
in Fig. 9. Figure 9 also shows that by increasing (3 from 0 to 0.3 and 0.5, the spurious oscillation is suppressed or 
eliminated. 
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Figure 1 0. One-dimensional Sod shock tube problem ; (a) and (b): comparison of exact and numerical solutions of density; (a) att — 2.0; (b) at 
t = 3.2, (b) shows how a shock, and a rarefaction (on the left hand side) are absorbed at the boundaries. The jumps at the two endpoints indicate the 
locations of the absorbing layers. 


5 Numerical Examples 

The present ETAU finite volume scheme is based on an unstructured grid topology and hence is flexible for com- 
plex geometries. With the dissipation model added to the basic scheme, as described in the previous section, the ETAU 
scheme covers a broad spectrum of flow from hypersonic to low subsonic, including aeroacoustic computations. Hence 
it is a time-accurate scheme of “all speed.” The robustness of the scheme is demonstrated in several numerical exam- 
ples below. The results are compared to either exact solutions, experimental data, or results from another numerical 
scheme. 


5.1 One-dimensional Riemann problem 

As a one-dimensional example, the Sod shock tube problem [25] is considered. This example also demonstrates 
how a shock and a rarefaction wave exit the one-dimensional domain without spurious reflection. For x G [0,8], 
initially, the two flow states are separated at x = 3: 


(p,u,p) 


(1,0,1), ifx < 3; 
(.125,0,0.1) ifx > 3 ' 


The physical boundary conditions/NRBCs at the left and right ghost cell centers, x = —0.01 and x = 8.01, are constant 
and respectively set at 

V L = (l,0,l) r , and V s = (0.125, 0,0.1) r . 

There are 400 uniform cells in [0, 8]. The ETAU scheme is employed with time step size At = 0.004. Figure 10a and b 
show respectively comparisons of computed and exact results before ( at t = 2.0) and after ( at t = 3.2) the rarefaction 
and the shock exit the boundary point x = 0 and x = 8. The shock and the rarefaction are all well captured. A van 
Albada limiter with a = 2 and (3 = 0 are employed. The effectiveness of the NRBCs is also demonstrated. The jumps 
near the two end points represent the absorbing layers, which are typical for the absorbing NRBC. 


5.2 Shock reflection and multiple shock-vortex interactions 

Consider Yee’s problem in a rectangular domain with 0 < x < 400, and 0 < j< 100 as shown in Fig. 11. A 
uniform triangulated grid of 40,000 cells is employed. A supersonic flow with a Mach number of 2.9 is given as the 
inflow boundary condition 

w 0 = 2.9, v 0 = 0, po = 1/1.4, po = l. 
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The boundary condition at the top is an inclined flow: 

u top = 2.6193, v top = -0.50632 ,p top = 1.5282, p top = 1.7000. 

The outflow boundary condition is the extrapolation NRBC and the bottom is a solid reflecting wall. Then, a steady 
oblique shock is formed with 29° inclination and reflected at the bottom wall. Figure 1 1 demonstrates the density 
contours and a comparison of the computed density and the exact solution along the centerline. Here a van Albada 
limiter with a = 2 and P = 0 (no dissipation model) are employed. 

This steady flow is then used as the initial background mean flow for a further computation of vortex- shock inter- 
action [26] . Such computation requires the scheme to handle both acoustic waves and strong shocks simultaneously. 
At t = 0, a strong Lamb vortex is placed at (22,60) (Fig. 12). With At = 0.05 and a = 0.5 in the van Albada lim- 
iter, 2640 time steps were computed. Figure 12 illustrates the multiple shock- vortex interactions at various times, 
and shows that the nonlinear acoustic waves are generated, and how they pass through the shocks and convect down- 
stream. Although there is no exact solution or experimental data, the side-by-side comparison in Fig. 12 shows that 
the numerical results by the present ETAU scheme and by the CE/SE scheme [3] are practically identical. 


5.3 Propagation of linear acoustic pulse and vortex/entropy waves 

Despite the capability of the present ETAU scheme in capturing discontinuities in flows of supersonic or hyper- 
sonic speed without suffering from the pathological behaviors, it is interesting to check its performance for flows at 
low Mach number with delicate acoustic waves. Here is an example illustrating the propagation of three basic types 
of weak, linear waves in a two-dimensional domain [26, 27]: linear acoustic pulse, vorticity and entropy waves. The 
computational domain in the x-y plane is a square with — 100 < x < 100, and —100 <y< 100. A uniform grid with 
40,000 uniform triangulated cells is used. Initially, a gaussian acoustic pulse is located at the center of the domain 
(0,0) and a weaker entropy /vorticity disturbance is located off the center at (67,0): 

p = 1 + 5 e -“i(-* 2 +3' 2 ) ; p = i + § e -“i(^ 2 +y 2 ) +0.l8e“ a2 K*- 67 ) 2+ A 

u = M + 0M8ye- a ^ x - 61 f +y2 l v=~0M8(x-67)e- a ^ x ~ 6 ^ 2+yl \ 

where y = 1.4, ai = ln2/9, and OC 2 = ln2/25. With a small amplitude factor, 8 = 0.001, the Euler equations are 
practically linearized. At all four boundaries, flow variables are the given M = 0.5 mean flow conditions: 


P = 



1, u = M = 0.5, v = 0. 
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vortex/shock interaction, time-accurate upwind, 
dt=0.05, alpha=0.5, dt between plots: 12, 








vortex/shocks interaction, CE/SE scheme, 
alpha=0.5, dt=0.05, dt between plots:12, 


Figure 12 . Vortex-shock interaction, isobars at different time steps, showing how non-linear acoustic waves are generated; top: by the current 
ETAU scheme ((3 = 0 ); bottom: by the CE/SE scheme [3], showing a side-by-side comparison. 
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Figure 13. Propagation of linear waves ; left: instantaneous density contours at time t = 63, right: comparison between exact and numerical 
solutions along the centerline y = 0 att = 63. When the entropy wave is exiting, notice that no spurious reflection is observed. 



Figure 1 4. Aeolian noise problem for a single cylinder; snapshot of isobars left: and isomachs right: at time step 250,000. 


These boundary conditions also play the role as an absorbing NRBC. Figure 13 shows the density contours at time 
t = 63, and its comparison to the exact solution along the v axis at t = 63. It is observed that the numerical solution 
agrees well with the exact solution and that no reflection is seen along the v axis. For this short term running, (3 = 0 is 
chosen. 


5.4 Aeolian noise of flow past single and twin cylinders 

This example is an aeroacoustics computation which demonstrates the time-accuracy of the scheme at low Mach 
number and for viscous flows. Aeolian noise is the noise generated by flow past a circular cylinder (or cylinders). 
First, consider a single circular cylinder with a diameter D = 1 .9 cm subject to a flow with a free stream Mach number, 
M = 0.2 [28]. The diameter, D , is used as the length scale. For convenience, the free stream speed of sound and density 
are used as scales for nondimensionalization. The Strouhal number (nondimensional frequency) is here conveniently 
defined based on D and the free stream speed of sound. The computational domain contains about 91,000 triangular 
cells, spanning between —2.5 < v < 10, and — 5 < y < 5. With a van Albada limiter (a = 0.5) and P = 0.0001, 80, 000 
time steps were computed to ensure the initial start-up transients have exited the domain. During the next 170,000 
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PSD at (-1,2), beta=0.0001 


PSD at (-1,2), beta=0.3 



Strouhal No. Strouhal No. 


Figure 15. PSD at a field point (—1,2) with different (3s; left: (3 = 0.0001; right: (3 = 0.3; showing that (3 has no influence on the PSD and 
Strouhal No. at the peak (PSD=— 4.36, St. =0.037). Due to the larger (3, the PSD curve on the right decays faster with increasing Strouhal No. 


Table 1 . Comparison of experimental and computed aeolian noise frequencies 


Mach No. 

Reynolds No. 

No. of cylinders 

Exp. Freq. 
(Strouhal No.) 

Comput. Freq. 
(Strouhal No.) 

0.2 

90,000 

single 

0.0369 

0.037 

0.0714 

15,800 

twin 

0.0146 

0.015 


time steps the pressure history at an observation point ( — 1 , 2) in the field are recorded. The point is located outside the 
vortex sheets to ensure that only acoustic wave fluctuations are recorded. Figure 14 shows the instantaneous isobars 
and Mach number contours at this final time step. Vortex sheets downstream of the cylinder are clearly observed. 
Then, fast Fourier transform (FFT) analysis is applied to the pressure history. The results are presented as a PSD 
(power spectral density) plot in Fig. 15. There are 17,000 sampling points (recorded every 10 time steps) in the 
pressure history file, which is just enough to cover 2 14 = 16, 384 points for FFT. In Fig. 15, the single peak of PSD 
corresponds to a Strouhal number, St = 0.037, which compares very well with the experimental data [28] (Table 1). 
Figure 15 also shows that the value of (3 has little influence on the Strouhal numbers at the curve peaks. However, the 
larger value of (3, e.g., 0.3, represents larger dissipation, the corresponding PSD curve is smoother and decays faster 
for St > 0.2. 

Similar computations are conducted for the aeolian noise problem with twin cylinders. The diameter of the 
cylinders is D = 0.955cm and the cylinders are placed vertically 3D apart (Fig. 16) and subjected to a freestream with 
a Mach number, M = 0.0714 [29]. The computational domain is about the same size as the single cylinder case, with 
a small buffer zone at the downstream outflow boundary. There are about 100,000 triangular cells in the domain. At 
a designated point (-1,2), the pressure history is recorded after 80,000 steps at every 20 steps for 660,000 additional 
time steps. There are 33,000 sample points in the pressure history file which slightly exceeds 2 15 = 32,768 for FFT. 
Table 1 demonstrates that for either single or twin cylinder, the computed frequencies (Strouhal numbers) agree well 
with the experimental data [28, 29]. Thus the present ETAU scheme is appropriate for viscous flow at very low Mach 
number. Figure 16 shows the instantaneous isobars, isomachs at time step 740,000 and the interaction between the 
two vortex streets from the twin cylinders. 
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Figure 16 . Snapshots of isobars left: and isomachs right: at time step 740,000 for the twin cylinder aeolian noise problem at free stream 
M = 0 . 0714 , (3 = 0 . 001 , showing vortex streets as well as the grid, the buffer zone is not shown. 


6 Concluding Remarks 

As a result of striking a balance between accuracy, efficiency, and affordable computer resources (CPU time 
and memory), an enhanced time accurate upwind (ETAU) finite volume method for unsteady Euler/Navier-Stokes 
equations is presented. The ETAU scheme is nominally second order accurate in space and time. The scheme adopts 
concepts from high order upwind schemes, such as the Cauchy-Kowalewski time stepping, accurate evaluation of 
surface fluxes, multidimensional limiter, etc. The Roe approximate Riemann solver is employed for its accuracy and 
efficiency. Also, employment of the unstructured grid enables flexibility in geometry. A built-in dissipation model in 
the ETAU scheme helps to overcome the usual pathological behaviors of the upwind schemes. The parameter (3 in the 
model is associated with the amount of artificial dissipation imposed. The computed flows are generally insensitive to 
the (3 values. Usually, P = 1CU 4 — 10 -3 is appropriate for most of the flows. Only when computing flows of different 
regimes (e.g. supersonic flow with high Mach number) does P need to be increased. 

With the dissipation model, the ETAU scheme is robust and viable for practical computations. It is a scheme of 
almost “all speed,” from Mach number 10 to the level as low as 10 -2 . It is time-accurate and works well with viscous 
or inviscid flows that are associated with either strong shocks or with delicate acoustic waves. The multidimensional 
limiter in the scheme also helps to establish a simple nonreflecting boundary condition. 
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